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Abstract 

We show that the thermodynamic hmit of a bistable phosphorylation-dephosphorylation cycle 
has a selection rule for the "more stable" macroscopic steady state. The analysis is akin to the 
Maxwell construction. Based on the chemical master equation approach, it is shown that, except 
at a critical point, bistability disappears in the stochastic model when fluctuation is sufficiently 
low but unneglectable. Onsager's Gaussian fluctuation theory applies to the unique macroscopic 
steady state. With initial state in the basin of attraction of the "less stable" steady state, the 
deterministic dynamics obtained by the Law of Mass Action is a metastable phenomenon. Stability 
and robustness in cell biology are emergent stochastic concepts. 
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The statistical physics of a hving cell requires a theory for open molecular system with 
chemical driving force and free energy dissipation Such system is capable of reaching a 
self-organizing state to which many biological functions are attributed. The state has been 
widely known as nonequilibrium steady- state (NESS) following M. Klein's concise terminol- 
ogy [2]. Two of the most exciting recent developments in statistical physics are concerned 
precisely with the NESS: The fluctuation theorem studies the novel entropy production 
characteristics of a NESS 31; and the one-dimensional exclusion process deals with highly 
non-trivial phase behavior 4|. 

However, the concept of NESS requires further refinements. This is the objective of this 
letter. From a statistical mechanics perspective, a NESS is a fluctuating, stochastically sta- 
tionary process. It has a stationary probability distribution as well as correlation functions 



For a wide class of physical and biological systems, this state is unique [6|. However, 
from a macroscopic perspective, an open, driven system can have multiple steady states. In 
fact, the dynamics can be even more complex that include oscillations, chaotic dynamics, 
and spatial-temporal chaos 

"Macroscopic" studies of living cell biochemistry are usually based on deterministic non- 
linear differential equations according to the Law of Mass Action 8| . Currently, it is generally 
accepted that a bistability in the deterministic dynamics corresponds to a bimodal proba- 
bility density function in the stochastic approach [o], [l^. With increasing size of a chemical 
reaction systems Uj, there is a separation of time scale: The transition rates between the 
two "macroscopic" states become infinitesimal ~ e~"^, where V is the systems volume and 
a is a positive constant. See [l^ for a detailed exposition. 

On the other hand, there is a well-developed, phenomenological fluctuation theory of 
NESS in statistical physics, pioneered by Onsager and Machlup, Lax, and Keizer, among 
many others [l^. One of the most important conclusions from this classic NESS fluctua- 
tion theory is a multivariate Gaussian fluctuations around a NESS. This result essentially 
conforms with Einstein's equilibrium fluctuation theory. 

In this letter, we shall use a concrete example to provide insights into this seeming paradox 
between the current view of NESS fluctuation, with multiple macroscopic steady states, 
and the classic Einstein-Onsager-Lax-Keizer (FOLK) Gaussian theory. Using the chemical 
master equation (CME) as the tool and a bistable system from current biochemical literature, 
we show that as ^ oo, a Maxwell-like construction is necessary. Such a construction 
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effectively singles out one unique state (or an attractor in the case of more complex dynamics) 
in the thermodynamics limit. It is around this state that the EOLK theory applies. Our 
analysis confirms the claim that the information necessary for the Maxwell-type construction 



is not present in the deterministic differential equation model of the system [13| ; one requires 
to build a mesoscopic, mechanistic model with stochasticity in order to gain the required 
information. Our conclusion is that, when dealing with biochemical reaction systems, one 
needs to differentiate the thermodynamic limit of a mesoscopic system and the differential 
equations based on the Law of Mass Action. The later follows the Kurtz's theorem 
the thermodynamics limit, however, has to augment the Maxwell construction, based on 
which the EOLK fluctuation theory applies. 

Several papers have addressed related issues in the past. We choose to revisit this impor- 
tant and fundamental issues in NESS due to the recent resurgent interests in the CME and 
,t. app.,catio„. to ceUula. bioche^str. In add.tion to Q. Ke.e. developed a Maxwell- 
type constructions for multiple nonequilibrium steady states [ig!]. While the present paper 
shares a similar idea, the previous approach was based on the diffusion approximation to a 
CME, an approach that can fail to represent correctly the mesoscopic steady-state , now 
known as Keizer's paradox icl|. 

The analysis performed here is for a particular example of a biochemical cycle, but what 
we show is general for nonlinear, driven chemical reaction systems with multistability. 

Phosphorylation- dephosphorylation cycle and the CME. Biochemical infor- 
mation processing inside cells uses a canonical reaction system called phosphorylation- 
dephosphorylation cycle (PdPC) flS^. We consider the Ferrell's kinetic model for PdPC 

18 1, which includes a positive feedback step, and its reversible extention first studied in 

19|: 



E + ATP + K* ^ E* + ADP + K*, 

a_i 

K + 2E*^K*, E* + P^E + Pi + P, (1) 

a_3 a-2 

in which E and E* are inactive and active form of a signaling protein. K and P are 
enzymes, kinase and phosphatase, that catalyze the phosphorylation and dephosphorylation 
respectively. K and K* are active and inactive forms of the kinase. The chemical reaction 
of ATP hydrolysis ATP ^ ADP + Pi provides the chemical driving force of the reaction. In 
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fact the free energy from the reaction is AG = fc^Tln {aia2[AT P] / {a-ia-2[AD P][Pi])} . In 
a cell, K and P, ATP, ADP, and Pi are all at constant concentration, and [E] + [E*] = etot- 
[z] denotes the concentration of the species z. 

The system in ([T]) exhibits bistability according to a deterministic analysis based on the 
Law of Mass Action [18]. It has been further shown that the bistability is distinctly a driven 
phenomenon that requires a sufficient large free ener gy d issipation 



its mesoscopic stochastic model in terms of the CME 



15| 



19|. Here we consider 



We shall denote ki = 0103 [ATP] /a_3, = a-ia^ADP]/a_^, /c2 = ci2[P] and k_2 = 
a^2[Pi][P]- Following the previous treatment flflQ, we assume the reversible binding 
K + 2E* ^ K* is rapid. The model thus is simplified into: E ^ E* with forward and 
backward rates = {kilKjx"^ + k^2){.^tot — 2^)5 R~{x) = {k2 + k_i[K]x'^)x, and x(t) = 

[E*]{t). The energy parameter from ATP hydrolysis 7 = exp{AG/kBT) = kik2/{k-ik-2)- 
7 = 1 is equivalent to a non-driven system which reaches a unique equilibrium steady state. 
In fact, the equilibrium probability distribution for the number of E* is binomial. 

The deterministic kinetic model based on the Law of Mass Action is 

^ = R+{x) - R- (x) = r(x; 6, e) (2) 

= k2 {9x^ [{etot -x) - ex] + [fi{etot ~ x) - x]] , 

in which the three parameters 6 = ki[K]/k2 represents the ratio of the activity of the kinase 
to that of the phosphatase; e = k^i/ki represents the ADP to ATP concentration ratio, 
and fi = k_2/k2 represents the strength of phosphorolysis. In a living cell, both fi and e are 
small; hence 7 = l/(/Ue) ^ 1. 

The fixed points of the Eq. ([3]) are the solution to R^{x) = R~{x). Their stability are 
determined by the -^{R^ix) — R~{x)). For some parameter ranges, the Eq. ([3]) exhibits 
saddle-node bifurcations 13], as shown in Fig. [1] One obtains the parameter region for 
the bistability from simultaneously solving r{z) = and = 0, which gives the boundary 
of the region of bistability, with a cusp, in {6, e) space (in terms of z as a parametric curve): 
Q ^ 2(1 + _ 3/i ^ ^ 2fielt - (/X + l)zetot _ ^ 

For the stochastic model in terms of the CME, one is interested in the number of E*, 
X, rather than its concentration. X takes non-negative integer values and is related to 
X = X/V where V is the system's volume. While X{t) is stochastic, its probability, P{X,t) 



satisfies the CME 
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Bifurcation diagram for steady states 
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FIG. 1: Bifurcation diagram of the steady states of reaction system in ([T]) and Eq. ([3]), 
function of the parameter ki. SN denotes saddle-node bifurcation. Other parameters used in the 
calculation: [K] = 1, etot = ^-i = 0.01, /c2 = 10, k-2 = 0.5. 



^ = VR^ (^) P(X - 1, t) + VR- (^) 
xP(X + l,t)-V [R+ ) + R- (f )] P{X,t). 
Its stationary solution gives the probability distribution in the NESS: 

P-(X) = P,(0)n ^""^ 



(4) 



Maxwell-type construction and stochastic bifurcation. When V is large, by the 
Euler-MacLaurin summation formula 



(5) 



where 



R-iy) 

etot ln(etot -x) -x\n 



+2\l ^ arctan 

u 



x{eex^ + 1) 
2 



-X I — 



arctan v 6'ex 



(6) 



We note that 



d(f){x) 
dx 



\og{R+{x)/R-{x)) 



In 



{etot - x){9x^ + /i) 



xidex'' 
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Hence the two stable fixed points of Eq. (j3j) correspond to tlie two minima of 4>{x), and the 
unstable fixed point corresponds to a maximum. In fact, for each steady state x*, 



<P'\x*) = 1 

X* 



R+jx) 

dlogx 



(7) 



x=x* 



which has the same sign as d{R — R^){x*) / dx. x* is stable if (j)' {x*) > 0, and unstable 
otherwise. Near a stable x* 

= <f>{x*) + (x - x*f + ■■■. (8) 

The Gaussian variance of P'^^{x) is which tends to zero when V tends to infinity. 

The square bracket term in ([7]) is called elasticity [20I due to its analogue to classical 
mechanics. Near the "more stable" stable fixed point of Eq. ([3]), for system with large V , a 



Gaussian, linear approximation is warranted. This is the classic theory of EOLK [12| . The 
key insight of this theory is the so-called the fluctuation-dissipation theorem for the NESS, a 
consequence of the Markovian Gaussian process. It provides a relationship among the linear 
relaxation kinetic matrix, the noise amplitude, and the covariance matrix of the Gaussian 



process. In a similar spirit. Berg et al. have put forward a linear noise approximation 20|. 



2l| has further illustrated that Gaussian characteristics is not necessarily only related to 



equilibrium fluctuations. Rather, it is determined by linear dynamics near a steady state 
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To our current discussion, the most important feature of Eq. ([5]) is that the function 
0(x) is independent of V , provided that V is sufficiently large. Therefore, even though 0(a:) 
exhibits bistability which corresponds closely to the Eq. ([3]), when ^ ^ 00, only one of 
the two stable fixed points is relevant in the thermodynamic limit, and it is the one with 
smaller 0(a;). A Maxwell-like construction, therefore, is necessary at the critical ki when 
<p{xl) = See Fig. [2 

It should be emphasized that bistability in the CME is really a nonequilibrium phe- 
nomenon. The metastability, however, can exist even when the white noise is "additive", 
or example, a diffusive particle restricted in a bistable potential with vanishing diffusivity 
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Discontinuity of stochastic entropy production and first-order phase transi- 
tion. A biochemical NESS is sustained by a continuous input of chemical energy which is 
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(b) Approximated distribution of CIVIE 
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(c) Exact distribution of CME 
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(d) Maxwell construction 
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FIG. 2: Stationary distribution of CME and Maxwell construction, (a) The shape of varies 
with the parameter ki; (b) Approximated stationary distribution of CME according to Eq. [5] 
varies with the volume with parameter value ki = 50; (c) Exact stationary distribution of CME 
according to Eq. U] varies with the volume with parameter value ki = 50; (d) The saddle-node(SN) 
bifurcation diagram for steady states and the Maxwell construction(MC). Other parameters used 
in the calculation: [K] = 1, etot = ^-i = 0.01, k2 = 10, A;_2 = 0.5. 



converted to dissipated heat: Entropy is produced in the process. The entropy production 
rate (epr) can be computed [lo| : 



AG 

fOD 




X 



+Jf f e-^'^^^'^dx ) . (9) 

Jx'^-e J 

The expression here is only valid when the reaction is sufficiently fast compared to diffu- 
sion, so that the reaction rate is only depend on time and not on the reaction coordinate. 

When V tends to infinity, epr/V tends to AGJf"* if (f){xl) < (p{x2), and tends to AG J. 
if 0(x2) < 0(3^1)- Since J^^ 7^ J|*, the epr/V is discontinuous at the critical situation when 

(j){xl) = 



ss 
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In classical equilibrium phase transition theory, a first-order phase transition has a dis- 
continuity in the first derivative of the free energy, and a second-order phase transition 
has a discontinuity in the second derivative. According to this classification, the present 
(nonequilibrium) phase transition can be considered as first order. 



Summary. A state of a biological cell, called a. functional cellular attractor [10], should 
be dynamically stable against various minor purturbations which are inevitable in living 
systems. Thus, it is often thought that "noise" added to the biological models only provides 
moderate refinements to the behaviors otherwise predicted by the classical, deterministic 
description. The present letter, however, shows something deeper: The relative stability and 
robustness of the phosphorylation-dephosphorylation module can not be properly inferred 
without an explicit consideration of the intrinsic noise in the model. In cellular biology, it is 
incorrect to model biological stability and robustness in terms of deterministic trajectories or 
sizes of basins of attractors from a deterministic model. Biological stability and robustness 
are stochastic concepts. Hence the presence of noise not only leads to corrections to the 
deterministic analysis but may give rise to emergent behaviors. 

The CME has now been recognized as a fundamental mathematical theory for meso- 
scopic chemical and biochemical reaction systems in a small, spatially homogeneous volume 

n 

|15l |. Its large volume limit recovers the Law of Mass Action kinetics [14]. However the 
deterministic differential equations, while define various attractors, provide no information 

n 

on the relative probabilities between them pjj. Furthermore, in the thermodynamic limit 
only one of the attractors will be dominant with probability 1. The Maxwell-type construc- 
tion, thus, enters the CME and becomes an integral part of a more complete theory. The 
biochemically interesting emergent dynamics from a CME, thus, is not the deterministic 
differential equations, but rather a stochastic jump process within a set^ discrete states 
defined by the deterministic attractors. This is distinctly a mesoscopic [24] driven system 
phenomenon: When the volume is too large, the time of escaping an attractor is practically 
infinite. Thus, the complex dynamics disappears. When there is no chemical driving force. 



i.e., . ^ 1, the n,ult.tabiMy d.appea. M- Nea. a given att.acto. wMch . a dete..i,n.tic 

fixed point, the EOLK phenomenological Gaussian fluctuation theory applies p^. Further- 
more, macrosccopic driven system in NESS can behave like an equilibrium system with a 
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(non- gradient) potential 



25] 
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